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A  CONTINUUM  MODEL  FOR  STREAMFLOW  SYNTHESIS 


1.  INTRODUCTION 

Streamflow  evolves  as  a  continuum,  and  is  normally  comprised  of  three 
components:  (1)  surface  runoff,  (2)  interflow,  and  (3)  baseflow.  These 
components  occur  concurrently,  although  their  relative  magnitudes  vary  with  time.  If 
we  consider  a  sudden  burst  of  rainfall,  then  surface  runoff  predominates  during  the 
rising  part  of  the  streamflow  hydrograph,  interflow  predominates  during  the  early  part 
of  its  recession,  and  baseflow  predominates  during  the  delayed  part  of  its  recession. 
The  mechanisms  and,  therefore,  the  governing  equations,  of  these  components  are 
different  but  are  influenced  by  dynamic  interactions  prevailing  between  them. 

Although  streamflow  synthesis  has  long  been  a  subject  of  scientific  inquiry, 
treatment  of  streamflow  as  a  continuum  taking  into  account  dynamic  interactions 
amongst  its  components  has  not  yet  been  fully  developed.  Most  of  the  approaches 
of  streamflow  synthesis  are  based  on  the  concepts  embodied  in  Horton’s  infiltration 
theory  of  runoff  (Horton,  1933).  According  to  this  theory,  rainfall  is  absorbed  for 
intensities  not  exceeding  infiltration  capacity,  while  for  excess  rainfall  there  is  a 
constant  rate  of  absorption  as  long  as  the  infiltration  capacity  is  unchanged.  Thus, 
infiltration  divides  rainfall  into  two  parts.  One  part  travels  over  the  surface  giving  rise 
to  surface  runoff,  and  the  other  part  infiltrates  into  the  ground  resulting  in 
replenishment  of  soil  moisture  and  recharge  of  groundwater,  and  eventually  in 
interflow  and  baseflow. 
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The  three  components  of  streamflow  have  been  treated  at  various  levels  of 
mathematical  sophistication  but  in  virtual  isolation  with  one  another.  Surface  flow 
has  been  studied  for  over  half  a  century  (Woolhiser,  1982),  and  as  a  result,  it  is 
understood  reasonably  well  (Hall,  1982).  Likewise,  baseflow  contribution  to 
streamflow  has  been  studied  for  nearly  30  years  and  it  too  is  understood  reasonably 
well  (Hall,  1982).  The  same,  of  course,  cannot  be  said  about  interflow.  This  is  not 
even  well  defined  and  is  least  understood.  Also,  least  understood  are  the  dynamic 
interactions  prevailing  amongst  these  components. 

Although  considerable  progress  has  been  made  in  mathematical  and 
numerical  treatment  of  the  boundary  value  problems  dealing  with  flows  over 
impervious  beds,  the  understanding  of  surface  flows  over  porous  beds  which 
dynamically  interact  with  subsurface  flow  is  quite  limited.  The  importance  of  this 
interaction  has  been  pointed  out  in  the  past  in  the  context  of  border  irrigation 
(Parlange,  1973),  and  in  the  study  of  flood  waves  in  ephemeral  streams  (Smith, 
1972).  These  studies,  however,  are  not  based  upon  a  coupled  set  of  equations 
pertaining  to  surface  and  subsurface  flow;  rather  the  attenuation  in  surface  flows  is 
included  by  considering  certain  infiltration  rate  with  time  lag.  The  dynamic  diffusion 
due  to  infiltration,  therefore,  remains  unaccounted  for. 

Freeze  (1972)  was  probably  the  first  to  develop  a  comprehensive  quantitative 
treatment  of  hillslope  hydrology  considering  explicit  interactions  between  near¬ 
surface  groundwater  flow,  surface  runoff  and  rainfall  intensity  patterns.  Rather 
limited  work  has  since  been  done  along  the  lines  of  Freeze.  Some  notable 
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examples  are  the  conceptual  model  of  Seven  and  Kirby  (1979),  the  model  of  Hillel 
and  Hornberger  (1979),  and  the  finite  element  model  of  Seven  (1977). 

The  most  recent  work  representing  a  major  step  forward  in  developing  an 
analytical  treatment  of  interdependent  surface  and  subsurface  hydrologic  processes 
is  by  Smith  and  Hebbert  (1983).  In  their  model,  the  hillslope  was  considered  to 
consist  of  two  soil  layers  with  the  lower  soil  capable  of  restricting  vertical  flow  at  the 
interface  creating  a  perched  aquifer  and  subsurface  stormflow.  Unsaturated  vertical 
flow  was  routed  by  a  kinematic  wave  method  and  linked  with  an  analytical  infiltration 
model.  Thus,  this  model  attempted  to  integrate  most  of  the  major  hydrologic 
response  mechanisms  presently  identified  as  contributing  to  the  hydrology  of  a 
simple  hillslope.  Other  hillslope  hydrological  models  (Gundy,  et  al.,  1985;  Stagnitti, 
et  al.,  1986)  and  surface  irrigation  models  (Walker  and  Humpherys,  1983;  Stagnitti, 
et  al.,  1986)  and  surface  irrigation  models  (Walker  and  Humpherys,  1983;  Ram,  et 
al.,  1983)  have  also  been  developed.  However,  none  of  these  models  developed  a 
method  to  compute  infiltration  rate  dynamically,  although  it  is  one  of  the  major 
factors  affecting  runoff  (or  advance  front)  and  surface  water  profile.  Most  of  the 
models  utilized  empirical  formulae  such  as  Kostiakov's  or  Green  and  Ampt's,  etc. 
Therefore,  the  prevailing  dynamic  process  between  surface  and  subsurface  flows 
remained  unaccounted  for. 

Toward  the  goal  of  eventually  accomplishing  a  continuum  model  for 
streamflow  synthesis,  three  related  areas  were  investigated:  (1)  subsurface 
unsaturated  flow,  (2)  flood  wave  propagation,  and  (3)  comparative  assessment  of 
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different  dynamic  wave  representations  of  shallow  water  wave  theory.  In  particular, 
we  focused  on  (a)  comparative  evaluation  of  the  kinematic-wave,  diffusion-wave, 
and  dynamic-wave  representations  of  the  shallow  water  wave  theory  ubiquitously 
applied  to  modeling  surface  runoff,  (b)  development  of  a  systems-  based  model  for 
infiltration,  and  (c)  modeling  movement  of  soil  moisture. 

2.  SURFACE  RUNOFF  MODELING 

The  surface  runoff  hydrology  was  investigated  along  three  lines:  (a) 
development  of  a  theory  of  errors  for  comparative  assessment  of  three  shallow 
water-wave  representations:  (i)  kinematic-wave,  (ii)  diffusion-wave,  and  (iii)  dynamic- 
wave;  (b)  development  of  discrete  linear  models  for  watershed  runoff;  and  (c) 
physically-based  Muskingum  methods  of  channel-flow  routing.  A  short  discussion  of 
each  is  in  order. 

2.1  Theory  of  Errors 

A  wide  range  of  problems  involving  free-surface  flows  can  be  modeled  using 
the  shallow  water-wave  (SWW)  theory.  The  SWW  theory  is  described  by  the  St. 
Venant  equations  or  their  approximations.  The  three  most  popular  representations 
of  the  SWW  theory  are  the  kinematic-wave  (KW),  diffusion-wave  (DW),  and  the 
dynamic-wave  (DYW)  approximations. 
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One  of  the  fundamental  questions  to  be  addressed  in  physically-based 
modeling  of  watershed  runoff  is  one  of  determining  the  appropriate  approximation  of 
the  shallow  water-wave  (SWW)  theory.  Of  the  different  approximations  of  the  SWW 
theory,  the  two  most  popular  approximations  are  the  kinematic-wave  (KW)  theory 
and  the  diffusion-wave  (DW)  theory.  How  accurate  are  these  approximations? 
Which  approximation  should  be  used  and  under  what  conditions?  What  is  the 
spatial  or  temporal  distribution  of  error  of  a  given  approximation?  What  is  the 
criterion  to  choose  between  these  approximations?  The  past  studies  have  dealt  with 
development  of  criteria  for  judging  the  adequacy  of  these  approximations.  However, 
these  criteria  are  point  values  and  do  not  relate  to  errors  resulting  from  use  of  these 
approximations.  Consequently,  the  error  in  space  and/or  time  is  not  known. 

The  larger  goal  of  this  study  was  to  develop  a  theory  of  errors  for  quantitative 
evaluation  of  the  adequacy  of  these  approximations,  and,  in  turn,  of  the  shallow- 
water-wave  theory.  However,  because  the  SWW  theory  consists  of  a  system  of 
nonlinear  partial  differential  equations  of  hyperbolic  type,  derivation  of  error 
equations  is  unattainable  at  this  stage.  Therefore,  some  realistic  simplifications 
were  made.  The  first  was  the  simplification  of  flows  being  time-independent. 

Steady  state  flows  are  ubiquitous  in  nature,  and  much  of  the  early  hydraulics  was 
based  on  this  assumption. 

Because  spatially  distributed  data  are  seldom  available,  it  was  assumed  that 
the  full  form  of  the  SWW  theory  or  the  dynamic-wave  representation  was  the  true 
representation  or  model,  and  was  capable  of  mimicking  the  behavior  of  the  real 
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world,  prototype  system,  and  the  kinematic-wave  and  diffusion-wave  approximations 
were  reasonable  approximations,  but  were  germane  to  conceptual  error.  The 
adequacy  of  these  approximations  is  well  documented  in  hydraulic  literature. 
However,  what  is  not  known  is  the  actual  error  and  its  distribution  in  time  and/or 
space,  as  well  as  its  relationship  to  flow  characteristics,  system  properties  and  initial 
and  boundary  conditions. 

The  theory  of  errors  can  serve  as  an  objective  criterion  for  judging  the 
adequacy  of  the  KW  and  DW  approximations  by  comparison  with  the  DYW 
approximation.  For  time-independent  flows,  the  theory  yields  error  as  a  function  of 
space  involving  infiltration,  and  boundary  conditions.  The  error  differential  equation 
is  ordinary  in  place  of  partial,  and  is  more  amenable  to  numerical  solution.  Even 
with  this  simplification,  analytical  solutions  are  not  possible  but  the  numerical 
solutions  are  much  simpler  and  easy  to  graph. 

Different  criteria  have  been  developed  to  evaluate  the  adequacy  of  the  KW 
and  DW  theories,  but  no  explicit  relations  either  in  time  or  in  space  between  these 
criteria  and  the  errors  resulting  from  these  approximations  have  yet  been  derived. 
Furthermore,  when  synthesizing  streamflow,  it  is  not  clear  if  the  kinematic-wave  and 
the  diffusion-wave  approximations  are  valid,  on  one  hand,  for  the  entire  hydrograph 
or  for  a  portion  thereof,  and  on  the  other  hand,  for  the  entire  channel  reach  or  for  a 
portion  thereof.  To  put  differently,  all  of  these  criteria  take  on  fixed  point-values  for 
a  given  rainfall-runoff  event.  This  study,  under  simplified  conditions,  derived  error 
equations  for  the  kinematic-wave  and  diffusion-wave  approximations  for  space- 


8 


independent  as  well  as  for  time-independent  flows.  These  equations  provided  a 
continuous  description  of  error  in  the  streamflow  hydrograph. 

With  these  considerations  in  mind,  errors  of  kinematic-wave  and  diffusion 
wave  approximations  were  derived  for  steady-state  channel  flows  subject  to  finite 
flow  at  the  upstream  end.  The  diffusion-wave  approximation  was  in  excellent 
agreement  with  the  dynamic  wave  representation  for  a  range  of  the  values  of  the 
Froude  number  and  the  kinematic-wave  number.  The  kinematic-wave  approximation 
was  also  in  good  agreement  with  the  dynamic  wave  representation,  but  for  a  limited 
range  of  the  values  of  the  Froude  number  and  the  kinematic-wave  number.  On  the 
other  hand,  the  approximate  diffusion-wave  analogy,  although  leading  to  analytical 
solutions,  was  not  accurate  and  should  not  be  employed. 

Under  two  different  initial  conditions  and  two  boundary  conditions,  solutions  of 
the  kinematic-wave  and  diffusion-wave  equations  were  derived  under  the 
simplification  that  the  flow  was  temporally  independent.  Thereafter,  error  equations 
for  the  KW  and  DW  theories  were  derived.  It  was  found  that  the  DW  theory  was 
quite  accurate  and  for  Froude  number  (Fo)  less  than  2  and  the  kinematic  wave 
number  (K)  greater  than  1 0,  the  DW  theory  would  be  an  accurate  representation  of 
the  SWW  theory.  Under  the  condition  where  there  was  no  downstream  control,  the 
KW  theory  was  an  accurate  representation  for  K  >  30,  K  F^^  >  5  .  The  KW  theory 
does  not  accommodate  a  downstream  control  and  hence,  as  expected,  did  not 
accurately  represent  the  SWW  theory  for  the  entire  channel  under  any  of  the  two 
boundary  conditions.  Details  of  this  work  are  contained  in  Singh,  Aravamuthan  and 
Joseph  (1994). 
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For  space-independent  flows,  a  dimensionless  parameter  was  defined, 
reflecting  the  effect  of  initial  depth  of  flow,  channel-bed  slope,  lateral  inflow, 
acceleration  due  to  gravity,  and  channel  roughness.  For  time-independent  flows,  the 
dimensionless  parameter  was  the  product  of  the  kinematic-wave  number  and  the 
square  of  the  Froude  number.  By  comparing  the  kinematic-wave  and  diffusion  wave 
solutions  with  the  dynamic-wave  solution,  error  equations  were  derived  in  terms  of 
the  aforementioned  dimensionless  parameters.  The  error  equations  for  space- 
independent  flows  turned  out  to  have  the  form  of  the  Riccati  equation.  The  work  is 
described  in  Singh,  Aravamuthan  and  Joseph  (1993). 

2.2  Watershed  Runoff 

Discrete  linear  models  were  developed  for  estimating  runoff  and  sediment 
discharge  hydrdgraphs  from  agricultural  watersheds.  A  regression  equation  was 
also  established  relating  runoff  rate  and  sediment  discharge.  Tested  on  five  small 
basins,  the  results  were  in  good  agreement  with  observations.  For  the  discrete 
linear  transfer  runoff  model,  the  values  of  the  integral  square  error  (ISE)  were 
generally  less  than  1%  for  all  calibration  events,  and  around  10%  with  the  average 
value  of  9.36%  for  all  verification  events.  For  the  discrete  linear  transfer  sediment 
model,  the  calibration  coefficient  of  determination  R  for  all  five  basins  was  more 
than  97%,  and  the  verification  R  was  more  than  91%  with  an  average  of  94.3%. 
Details  of  this  work  are  described  in  Wang  et  al.  (1991). 
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2.3.  Physically  -  Based  Muskingum  Method 

Flow  routing  in  channels  was  investigated  using  the  kinematic  wave  theory  as  well 
as  the  Muskingum  method.  For  the  latter  method  parameters  were  derived  from  the 
St.  Venant  equations.  Preliminary  testing  showed  that  this  method  of  parameter 
estimation  made  the  Muskingum  method  more  accurate  than  any  of  the 
conventional  methods.  The  kinematic  wave  method  was  investigated  for  perennial 
as  well  as  ephemeral  streams.  This  method  can  be  extended  to  include  flood  wave 
propagation  due  to  dam  rupture.  This  work  is  more  fully  described  in  Wang  and 
Singh  (1992). 

3.  SUBSURFACE  FLOW 

Modeling  of  flow  of  water  in  the  unsaturated  zone  is  far  from  complete, 
especially  at  the  field  scale.  Two  lines  of  inquiry  were,  therefore,  launched.  First,  a 
systems  approach  was  developed  to  model  infiltration  and  soil  moisture,  which  holds 
promise  for  unifying  different  infiltration  models  reported  in  the  literature.  This 
approach  can  also  relate  parameters  of  one  infiltration  to  another.  The  second  type 
of  approach  pertained  to  application  of  the  kinematic  wave  theory.  This  approach 
has  the  advantage  of  coupling  plant  root  extraction  and  evapotranspiration  with  soil 
water  dynamics. 

3.1  Infiltration  Modeling 

A  general  infiltration  model  was  derived  using  systems  approach.  The 
models  of  Horton,  Kostiakov,  Overton,  Green  and  Ampt,  and  Philip  are  some  of  the 
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example  models  which  are  shown  as  special  cases  of  the  general  model.  An 
equivalence  between  the  Green-Ampt  model  and  the  Philip  two-term  model  was 
shown.  The  general  model  also  provides  a  solution  for  the  Holtan  model  expressing 
infiltration  as  a  function  of  time.  This  solution  of  the  Holtan  model  does  not  appear 
to  have  been  reported  in  the  literature.  A  first-order  analysis  was  performed  to 
quantify  the  uncertainty  involved  with  the  generalized  model.  The  general  infiltration 
model  contains  five  parameters.  Two  of  the  parameters  are  physically  based  and 
can  therefore  be  estimated  from  the  knowledge  of  soil  properties,  antecedent  soil 
moisture  conditions,  and  infiltration  measurements.  The  remaining  three  parameters 
can  be  determined  using  the  least  squares  method.  The  model  was  verified  using 
ten  observed  infiltration  data  sets.  Agreement  between  observed  and  computed 
infiltration  was  quite  good.  This  work  is  more  fully  described  in  Singh  and  Yu 
(1990). 

3.2  Movement  of  Soil  Moisture 

The  unsaturated  subsurface  flow  serves  as  a  link  between  surface  runoff  and 
groundwater  runoff,  and,  therefore,  occupies  the  central  position  in  the  streamflow 
continuum.  The  unsaturated  flow  region  is  the  upper  soil  matrix  which  also  is  the 
principal  source  of  interflow.  For  surface  flow,  infiltration  is  the  major  sink,  and  for 
groundwater  flow,  soil  moisture  percolation  is  the  principal  source  or  recharge. 

Much  of  the  mathematical  treatment  of  unsaturated  flow,  reported  to  date,  is  based 
on  the  Fokker-Planck  equation  or  Richards  equation.  Both  these  equations  are  of 
the  diffusion  type,  and  do  not  lend  themselves  to  analytical  solutions,  except  for 
overly  simplified  cases. 
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If  one  ignores  the  effect  of  diffusion  and  models  soil  moisture  movement 
using  the  kinematics  wave  theory,  then,  under  certain  simplifying  but  realistic 
conditions,  it  is  possible  to  derive  analytical  solutions.  This  premise  was  pursued  in 
this  project.  Currently  available  one-dimensional  kinematic-wave  models  assume 
absence  of  sink  terms.  In  other  words,  once  the  water  gets  infiltrated,  it  is  either 
retained  by  the  soil  or  moves  downward  to  recharge  the  groundwater.  This 
assumption  is  not  tenable,  especially  in  agricultural  or  forest  watersheds.  In  this 
project,  an  effort  was  made  to  include  a  sink  term  in  modeling  of  soil  moisture.  This 
sink  term  may  represent  removal  of  soil  moisture  by  vegetation  through  the  process 
of  evapotraspiration. 

Recognizing  the  difficulties  of  modeling  unsaturated  flow  using  the  Richards 
equation,  a  novel  approach  was  developed  in  this  project.  This  approach  was 
based  on  the  kinematic  wave  theory  in  which  a  unique  relation  is  hypothesized 
between  the  flux  and  the  flow  concentration  or  the  hydraulic  head. 

The  fundamental  assumption  underlying  this  theory  is  that  the  moisture 
movement  is  gravity-dominated  and  the  hydraulic  conductivity-soil  moisture 
relationship  is  single-valued,  i.e.,  it  does  not  experience  any  hysteretic  effects. 
Although  these  assumptions  are  not  strictly  valid,  they  do  provide  a  reasonably 
accurate  approximation.  Another  advantage  of  the  theory  is  its  simplicity  and  that  it 
allows  analytical  solutions  under  simplified  conditions.  Under  more  complex 
conditions,  numerical  solutions  are  easily  derived. 
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This  unique  relation,  when  coupled  with  the  continuity  equation  expressing 
the  conservation  of  mass,  gives  rise  to  a  first  order,  nonlinear  hyperbolic  equation. 
Under  simplified  initial  and  boundary  conditions,  analytical  solutions  of  this 
kinematic-wave  equation  are  tractable.  Following  this  tract,  the  soil  moisture 
movement  was  modeled  with  the  use  of  the  kinematic-wave  theory.  However, 
attention  is  to  be  focused  on  certain  aspects  of  the  kinematic-wave  modeling  that 
are  not  apparent  at  the  first  glance.  First,  because  the  time  history  of  the  moisture 
front  distinguishing  between  wet  and  dry  soil  is  unknown,  the  kinematic-wave 
formulation  of  soil  moisture  movement  results  in  a  free-boundary  problem.  Second, 
natural  watersheds  have  vegetation  either  seasonally  or  throughout  the  year. 
Inclusion  of  vegetation  in  modeling  of  soil-moisture  movement  gives  rise  to  an 
additional  free  boundary,  further  complicating  the  model. 

With  these  considerations  in  mind,  a  kinematic-wave  model  was  developed 
for  sihiulating  the  movement  of  soil  moisture  in  unsaturated  soils  with  plants.  The 
model  involved  three  free  boundaries.  Analytical  solutions  were  derived  when  the 
plant-root  extraction  of  moisture  was  at  a  constant  rate,  and  the  upstream  boundary 
condition  was  time  independent.  If  these  assumptions  are  waived,  then  numerical 
solution  is  the  only  resort. 

With  the  use  of  this  theory,  a  comprehensive  analytical  treatment  has  been 
developed  for  soil  moisture  movement  with  plant-root  extraction.  The  treatment 
involves  free  boundaries  and  to  our  knowledge  this  has  not  been  dealt  with  in  the 
literature  so  far.  This  work  is  described  in  Singh  and  Joseph  (1994). 
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